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We derive a correct first-order perturbation theory in electromagnetism for cases where an interface between 
two anisotropic dielectric materials is slightly shifted. Most previous perturbative methods give incorrect re- 
sults for this case, even to lowest order, because of the complicated discontinuous boundary conditions on the 
electric field at such an interface. Our final expression is simply a surface integral, over the material interface, 
of the continuous field components from the unperturbed structure. The derivation is based on a "localized" 
coordinate-transformation technique, which avoids both the problem of field discontinuities and the challenge 
of constructing an explicit coordinate transformation by taking a limit in which a coordinate perturbation is 
infinitesimally localized around the boundary. Not only is our result potentially useful in evaluating boundary 
perturbations, e.g. from fabrication imperfections, in highly anisotropic media such as many metamaterials, but 
it also has a direct application in numerical electromagnetism. In particular, we show how it leads to a sub-pixel 
smoothing scheme to ameliorate staircasing effects in discretized simulations of anisotropic media, in such a 
way as to greatly reduce the numerical errors compared to other proposed smoothing schemes. 



I. INTRODUCTION 



In this paper, we present a technique to apply perturba- 
tive techniques to Maxwell's equations with anisotropic ma- 
terials, in particular for the case where the position of an 
interface between two such materials is perturbed, general- 
izing an earlier result for isotropic materials (1]. In this 
case, the discontinuities of the fields at the interface cause 
many standard perturbative methods to fail, which is unfor- 
tunate because such methods are very useful for many prob- 
lems in electromagnetism where one wishes to study the ef- 
fect of small deviations from a given structure — not only do 
perturbative methods allow one to apply the computational 
efficiency of idealized problems to more realistic situations, 
but they may also offer greater analytical insight than brute- 
force numerical approaches. The corrected solution described 
in this paper should aid the study of interface perturbations, 
from surface roughness to fiber birefringence, in the context 
of anisotropic materials. Such materials have become in- 
creasingly important thanks to the discovery of "metamate- 
rials," subwavelength composite structures that simulate ho- 
mogeneous media with unusual properties such as negative 
refractive indices [2], and which may be strongly anisotropic 
in certain applications — for example, those involving spher- 
ical or cylindrical geometries |3] such as recent proposals 
for "invisibility" cloaks J4). Furthermore, we have recently 
shown that interface-perturbation analyses benefit even purely 
brute-force computations, because they enable the design of 
sub-pixel smoothing techniques that greatly increase the ac- 
curacy (and may even increase the order of convergence) of 
discretized methods [5], which are normally degraded by dis- 
continuous interfaces |6][7). Here, we show that our corrected 
perturbation analysis provides similar benefits for modeling 
anisotropic materials, where it yields a second-order accurate 
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smoothing technique (correcting a previous heuristic proposal 

There have been several previous approaches to rigor- 
ous treatment of interface perturbations in electromagnetism, 
where classic approaches for small Ae perturbations fail be- 
cause of the field discontinuities 0] [8] E3 ED El- One 
approach that was applied successfully to boundaries between 
isotropic materials is essentially to guess the correct form of 
the perturbation integral and then to prove a posteriori that it is 
correct [ Q. For isotropic materials, where there is some guid- 
ance from effective-medium heuristics [13], this was practi- 
cal, but the correct answer (below) appears to be much more 
difficult to guess for anisotropic materials. Another approach, 
which generalizes to the more difficult case of small surface 
"bumps" that are not locally flat, was to express the prob- 
lem in terms of finding the polarizability of the perturbation 
and then connecting it back to the perturbation integral via 
the method of images fl2l . For a locally flat perturbation be- 
tween isotropic materials, this process can be carried out an- 
alytically to reproduce the previous result from Ref. JT], but 
it becomes rather complicated for anisotropic media. Third, 
one can transform the problem into a statement about the co- 
ordinate system to avoid problems of shifting field disconti- 
nuities, by finding a coordinate transformation that expresses 
the interface shift iflQlfTTI . This approach, while powerful, has 
two shortcomings: first, finding an explicit coordinate trans- 
formation may be difficult for a complicated interface pertur- 
bation; and second, the resulting perturbation integrals are ex- 
pressed in terms of the fields everywhere in space, not just at 
the boundaries. Intuitively, one expects that the effect of the 
perturbation should depend only on the field at the boundaries, 
as was found explicitly for the isotropic case (HE]- ^ n m i s 
paper, we derive precisely such an expression for the case of 
interfaces between anisotropic materials, by developing a gen- 
eral new analytical technique for interface perturbations: we 
express the perturbation as a coordinate transformation, but 
using a coordinate transform localized around the perturbed 
interface, and take a limit in which this localization becomes 
narrower and narrower so that the choice of transform disap- 



2 



pears from the final result. 

In the following sections, we first formulate the problem 
of the effect of an interface perturbation more precisely, re- 
late our formulation to other possibilities, and summarize our 
final result in the form of eq. ([3J. We then derive quite gener- 
ally how to formulate the problem of interface perturbations 
in terms of a localized coordinate transformation, and show 
how this allows us to express the perturbation-theory integral 
as a sum of contributions around individual points on the in- 
terface. Next, we apply this framework to the specific prob- 
lem of a boundary between two anisotropic dielectric materi- 
als, and derive our final result. As a check, our perturbation 
theory is then validated against brute-force computations for 
a simple numerical example. Finally, we discuss the appli- 
cation of our new perturbation result to sub-pixel smoothing 
of discretized numerical methods, and show that we obtain a 
smoothing technique that leads to much more accurate results 
at a given spatial resolution. In the appendix, we provide a 
compact derivation and generalization of a useful result lfl4l 
relating coordinate transformations to changes in e and /.i. 



II. PROBLEM FORMULATION 

There are many ways to formulate perturbation techniques 
in electromagnetism. One common formulation, analogous to 
"time-independent perturbation theory" in quantum mechan- 
ics fl5l . is to express Maxwell's equations as a generalized 
Hermitian eigenproblem V x V x E = u 2 eE in the fre- 
quency u> and electric field E (or equivalent formulations in 
terms of the magnetic field H) [ 16 1, and then to consider the 
first-order change Aui in the frequency from a small change 
Ae in the dielectric function e(x) (assumed real and positive), 
which turns out to be ifTBI : 



Alu _ f E* ■ AeE d 3 x 
77 ~ ~ 2 J E* • eE d 3 x 



0(Ae 2 ), 



(1) 



where E and u> are the electric field and eigenfrequency of the 
unperturbed structure e, respectively, and * denotes complex 
conjugation. The key part of this expression is the numerator 
of the right-hand side, which is what expresses the effect of 
the perturbation, and this same numerator appears in a nearly 
identical form for many different perturbation techniques. For 
example, one obtains a similar expression in: finding the per- 
turbation A/3 in the propagation constant f3 of a waveguide 
mode 021; the coupling coefficient (~ J E* • AeE') between 
two modes E and E' in coupled-wave theory lfT8l [T9l l20l : 
or the scattering current J ~ AeE (and the scattered power 
~ J 3* ■ E) in the "volume-current" method (equivalent to the 
first Born approximation) lfl"2l |2~T1 l22l [23 1 . Eq. ([T]i also corre- 
sponds to an exact result for the derivative of uj with respect to 
any parameter p of e, since if we write Ae = J| Ap+ 0(Ap 2 ) 
we can divide both sides by Ap and take the limit Ap — > 0; 
this result is equivalent to the Hellman-Feynman theorem of 
quantum mechanics lfTl [T5ll . In cases where the unperturbed e 
is not real, corresponding to absorption or gain, or when one 
is considering "leaky modes," the eigenproblem typically be- 




Figure 1: Schematic of an interface perturbation: the interface be- 
tween two materials e a and e b (possibly anisotropic) is shifted by 
some small position-dependent displacement h. 



comes complex-symmetric rather than Hermitian and one ob- 
tains a similar formula but without the complex conjugation 
E4ll . Therefore, any modification to the form of this numera- 
tor for the frequency-perturbation theory immediately leads to 
corresponding modified formulas in many other perturbative 
techniques, and it is sufficient for our purposes to consider 
frequency-perturbation theory only. 

As we showed in Ref. (TJ, eq. ([T]l is not valid when Ae is 
due to a small change in the position of a boundary between 
two dielectric materials (except in the limit of low dielectric 
contrast), but a simple correction is possible. In particular, let 
us consider situations like the one shown in Fig.[T] where the 
dielectric boundary between two materials e a and e b is shifted 
by some small displacement h (which may be a function of 
position). Directly applying eq. ([IJ, with Ae = ±(e a — e b ) 
in the regions where the material has changed, gives an incor- 
rect result, and in particular Auj/h (which should ideally go to 
the exact derivative dio / dh) is incorrect even for h — ► 0. The 
problem turns out to be not so much that Ae is not small, but 
rather that E is discontinuous at the boundary, and the stan- 
dard method in the limit h — > leads to an ill-defined surface 
integral of E over the interfaces. For isotropic materials, cor- 
responding to scalar e a,b , the correct numerator turns out to 
be, instead, the following surface integral over the boundary 
ID: 



E* ■ AeE fx 



En 



\Dj 



h • dA, (2) 



where E|| and D± are the (continuous) components of E and 
D = eE parallel and perpendicular to the boundary, respec- 
tively, dA points towards e b , and h is the displacement of the 
interface from e° towards e b . 

In this paper, we will generalize eq. Q to handle the case 
where the two materials are anisotropic, corresponding to 
arbitrary 3x3 tensors e a and e b (assumed Hermitian and 
positive-definite to obtain a well-behaved Hermitian eigen- 
problem). In the generalized case, it is convenient to define 
a local coordinate frame (%x , x%, x§) at each point on the sur- 
face, where the x\ direction is orthogonal to the surface and 
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the (x2,X3) directions are parallel. We also define a contin- 
uous field "vector" F = (Di,E 2 , E 3 ) so that F\ = Z?j_ and 
F 2 .3 = Em. As derived below, the resulting numerator of 
eq. ([TJ, generalizing eq. is: 



F* • [r (e a ) - r (e b )] Fh dA, 



(3) 



where r(e) is the 3 x 3 matrix: 



r(e) 



£11 
£21 
Ell 



£22 
£32 



£11 

£2l£l2 



£11 

£3l£l2 




(4) 



which reduces to eq. (|2| when e is a scalar multiple e of the 
identity matrix. (Our assumption that e is positive-definite 
guarantees that en > 0.) 

We should note an important restriction: eq. (j2j) and eq. (|3]l 
require that the radius of curvature of the interface be much 
larger than h = |h|, except possibly on a set of measure zero 
(such as at isolated corners or edges). Otherwise, more com- 
plicated methods must be employed fl2l . For example one 
cannot apply the above equations to the case of a hemispher- 
ical "bump" of radius h on the unperturbed surface, in which 
case the lowest order perturbation is Alj ~ 0(h 3 ) and re- 
quires a small numerical computation of the polarizability of 
the hemisphere IfTZll . 



III. LOCAL COORDINATE PERTURBATIONS 

The difficulty with applying the standard perturbation- 
theory result ([T]) to a boundary perturbation is that, instead 
of a small Ae with fixed boundary conditions on the fields (to 
lowest order), we have a large Ae over a small region in which 
the field boundary discontinuities have shifted. However, we 
can transform one problem into the other: we construct a co- 
ordinate transformation that maps the new boundary location 
back onto the old boundary, so that in the new coordinates 
the boundary conditions are unaltered while there is a small 
change in the differential operators due to the coordinate shift. 
In expressing the problem in this fashion, we will present two 
key techniques. First, we employ a result from lfl4l . general- 
ized in the Appendix to anisotropic materials, that expresses 
an arbitrary coordinate transform as a change Ae and A/j, in 
the permittivity and permeability tensors, which allows us to 
directly apply eq. ([TJ. Second, unlike Refs. iflOl [TD . we do 
not wish to explicirty construct any coordinate transformation, 
since this may become very complicated for an arbitrary per- 
turbation in an arbitrary-shaped boundary. Instead, we ex- 
press the boundary shift in terms of a local coordinate trans- 
form, that only "nudges" the coordinates near the perturbed 
boundary, and in the limit where the region of this coordi- 
nate perturbation becomes arbitrarily small we will recover 
the coordinate-independent surface integrals Q and (J3J. 



A. Coordinate perturbations 

Suppose that in a certain coordinate system x we have elec- 
tric field E(x,t), magnetic field H(x, t), dielectric tensor 
e(x), and relative magnetic permeability tensor ju(x), satis- 
fying the Euclidean Maxwell's equations. Now, we transform 
to some new coordinates x'(x), with a 3 x 3 Jacobian matrix 

Cf defined by — i ^ L . In the new coordinates, the fields 
can still be written as the solution of the Euclidean Maxwell's 
equations if the following transformations are made in addi- 
tion to the change of coordinates: 



E' = (J 7 



h' = (j T y 



E, 



H 



e = 



Je J 1 
detj 



V = 



JuJ 1 

detj 



(5) 



(6) 



(7) 



(8) 



where J" denotes the transpose. This result is derived in the 
Appendix, generalized from the result for scalar e and fi from 
Ref. fl4l . 

Now, suppose the coordinate change is "small," meaning 
that S = 1 + AJ~, where the eigenvalues of AjT(x) are ev- 
erywhere O(S) for some small parameter 6. Then Ae(x') = 
e'(x') - e[x(x')] = 0(6) and similarly Afj, = O(S). There- 
fore, the solutions of Maxwell's equations will be nearly those 
of e and fi merely translated to the new coordinate locations, 
and the difference due to Ae and Afi can be accounted for, 
to 0(<5 2 ), by first-order perturbation theory. That is, general- 
izing eq. ([TJ to the case of anisotropic media with both e and 
fi, one finds by elementary perturbation theory for the gener- 
alized eigenproblem: 



Alu 



f [E* ■ Ae ■ Eq + 

/[E^-e-Eo + H* 
/[E^-Ae-Eo + HS 



Afj, ■ H ] d 3 x' 
[A ■ H ] d 3 x' 
A/x • H ] d 3 x' 



2 JE * ■e-E Q d s x' 



0(S 2 ) 



(9) 



where the "0" subscripts denote the solution for the unper- 
turbed system, given by e[x(x')] and /x[x(x')], i.e. e and /j, 
simply translated into the x' coordinates without transforming 
by the Jacobian factors. 



B. Interface-localized coordinate transforms 



Suppose that we have an unperturbed interface between two 
materials e a and e b that forms a surface So (i-e-, the points 
x G So), and we perturb it to a new interface S by a small 
perpendicular shift h(x) as depicted schematically in Fig. [T] 
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In order to investigate this boundary shift, we will perform a 
coordinate transform x'(x) that shifts S to S' — Sq. That is, 
in our new coordinates, the interface has not been perturbed, 
but the materials have changed by the Jacobian factors as de- 
scribed in the previous section. Moreover, we will construct 
our coordinate transform so that it is localized to the interface, 
i.e. so that x' = x far from Sq. In particular, we write: 



then write 



x' = x - h(x)L(x) 



(10) 



where £(x) G [0,1] is some differentiable localized function, 
equal to unity on the interface [L(S) = 1] and identically zero 
outside some small radius-i?/2 neighborhood of the interface 
(the support of L lies within this neighborhood), chosen so 
that |Vi| = 0(1/ R). Eq. ([To} is constructed so that x 6 S 
implies x' e 
mapped to So 



S' = Sq, causing the new interface S to be 
as desired. Thus, h(x) for x 6 S must be 
the perpendicular displacement from So to S. For x ^ S, 
h(x) should be some differentiable, slowly-varying function 
(except possibly at isolated surface kinks and discontinuities). 
The precise functions L and h will turn out to be irrelevant to 
our final answer ([3]), so we need not construct them explicitly. 

We will take |h| = h <C 1 to be the small parameter of our 
perturbation theory, and will concern ourselves with obtain- 
ing the correct first-order Acu in the limit h — > 0. We will also 
eventually take the limit R — ► 0, but will still require h ^ R 
in order to ensure, as shall become apparent below, that the Ja- 
cobian factor of the coordinate transformation remains close 
to unity. (That is, we let h go to zero faster than R.) Finally, in 
order to have h(x) be sufficiently slowly varying that we can 
neglect its derivatives compared to the derivatives of L(x), be- 
low, it will be important to require that the radius of curvature 
of Sq and S be much larger than h, except possibly at isolated 
points; otherwise, more complicated perturbative methods are 
required fl2l . 



C. Point-localized coordinate transforms 

The coordinate transformation ( p~0] > representing our bound- 
ary perturbation is localized around the perturbed interface, 
but is convenient to go one step further: we will represent 
the coordinate transform as a summation of coordinate trans- 
formations localized around individual points on the interface, 
by exploiting the concept of a partition of unity from topology 
ll25l . reviewed below. 

Consider the support of the function L(x) from above. 
This support is covered by the open set of spherical radius-i? 
neighborhoods of every point on the surface, and that cover- 
ing must admit a locally finite subcovering {J/'"* 1 }; that is, 
a subset of neighborhoods {LA a )} such that every point on 
the surface intersects finitely many neighborhoods U^ a \ and 
the union of the covers the support of L. There must 
also exist a partition of unity {0( Q )}: a set of differentiable 
functions (/>(") (x) € [0,1] with support C U^ a \ such that 
$ ( x ) = 1 everywhere in the support of L. We can 



L(x) = 



E^ (Q) (x) 



L(x) = 5>(«)(x), 



(11) 



where each if( Q )(x) = 0( Q '(x) L(x) e [0, 1] is a differen- 
tiable function localized to a small radius-i? neighborhood 
t/( Q ) of a single point on the interface. The Jacobian J~ of 
the coordinate transformation ( fT0| ) can then be written in the 
form: 



where 



(«) 



dx,< 



fci(x) i^(x) 



(12) 



(13) 



has support C {JW. 

The key advantage of this construction arises if we look 
at Ae = e' — e from eq. Assuming Aj7" is small 

and we are computing Ae to first-order, then we can write 
Ae = J2 a Ae( Q ) as a sum of contributions from each AS^ a ' 
individually, and similarly for Afi. Therefore, when comput- 
ing the first-order perturbation Ao; from eq. Q, we can write 
Auj = ~J2 a Au;( a ) as a sum of contributions Au'°' analyzed 
in each point neighborhood separately. This removes the need 
to deal with the complex shape of the entire boundary at once, 
and is the procedure that we adopt in the following section. 



IV. PERTURBATION THEORY DERIVATION 

In the previous section, we established several important 
preliminary results that allow us to express a boundary per- 
turbation, via coordinate transformation, as a sum of local- 
ized material perturbations Ae^ and Afi^ around indi- 
vidual points of the boundary. We will now explicitly eval- 
uate those contributions, taking the limit as the perturbation 
h — ► and the coordinate distortion radius R — > to obtain 
our coordinate-independent final result, eq. 

We therefore restrict our attention to a single neighborhood 
[A Q ) and the contribution from the corresponding term 
in the coordinate transformation. In this small neighborhood 
of radius R, we can take h(x) w to be a constant to 
lowest order in R. In this case, the interface is locally flat, and 
we can choose a local coordinate frame {x\, X2, Xs) so that 
x\ is the direction perpendicular to the interface at x\ = 0, 
with x\ < corresponding to e a and x\ > corresponding 
to e b , as shown in Fig. [2] In this coordinate frame h( Q ) = 
(h( a \ 0, 0), the Jacobian contribution AS^ simplifies to: 



AJ-, a) = -6 iX h^K\ a) + 0(h^R), 



(14) 



where Sn is the Rronecker delta and denotes 
dK^ /dxj. Since R > |h| by assumption and K ^ € [0, 1] 
is a smooth localized function with support of radius R, 
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Figure 2: Schematic of an interface perturbation as in Fig. |T| magni- 
fying a small portion of the interface where the surface is locally flat. 
A local coordinate frame (xi, X2,xs) is chosen so that x\ is perpen- 
dicular to the surface, and so that x\ = denotes the location of the 
perturbed surface S (shifted perpendicularly by h from the original 
surface So). 



can be constructed so that h^K^ — 0(h/R), i.e so that 
the derivatives are small. This will make J" close to unity and 
allow us to use the perturbation equation (|9j. 

We must now construct Ae^ to first order. Since J' = 
1 + X; q AJ" (q) , we obtain: 

1 + V h^K[ a) + 0(h 2 ) + 0(hR). (15) 



detj 



Combined with eq. ( 13 1, we can now evaluate eq. ^ for e', to 
lowest order, to obtain Ae = J2 a Ae(a) + °( h2 ) + 0(hR), 
with 



As 



(a) 



(a) 



y 

k 



This will contribute to |9]) via the integral: 



j(a) = f E * . Ag (a) . Ed 3 X; 

Ju(") 



(16) 



(17) 



where we have dropped the "0" subscript from the unper- 
turbed field E for simplicity. In order to simplify this in- 
tegral, we will write E = {E 1 ,E 2 ,E 3 ) in terms of F = 
(Di, E 2 , E 3 ), since F is continuous whereas E\ is not. Solv- 
ing for Ei inD = e-E yields E\ = j^{Di—e\ 2 E 2 —£i 3 E 3 ), 
and thus E = !F F where 



1 

EH 



. £12 

£11 

1 



£13 
Ell 



(18) 



Because F is continuous, we can write F(x) = + O(R), 
where the O(R) term is a higher-order contribution to I^ a ' 
that can be dropped and the F^ Q ^ is a constant that can be 
pulled out of the integral. Therefore, we are left with 



(7(«> 



j(a) =p (o)*. / T [ ■ Ae i::i ■ T.Tx ■ F !< " + 0(hR), 

(19) 

where H is the conjugate-transpose. This integral now sim- 
plifies a great deal, because the only non-constant terms are 



from the and the step-function 6(a;i) dependence of 

e(x). In particular, the integrals over the K^ a) and K { 3 a) terms 
vanish, because along the x 2 and x 3 directions, respectively, 
they are integrals of the derivatives of a function K ( a > that 
vanishes at the endpoints. We are left with the terms, 
which yield the integrand: 



T< ■ I £ 22 £23 I ■FK[ a) h^ = T(e)K[ a) hM 

£32 £33 / 

= {r(e a ) + [r(s b ) - r(e a )] Q(xi)} K<f ] h^ (20) 

where the product of the three matrices gives precisely the 
matrix r(e) defined in eq. Q, u sing the assumption that e is 
Hermitian (gl = e). When eq. ( |20| ) is integrated by parts in 
the x\ direction, we obtain the integral of K^ a > multiplied by a 
delta function 8{xx) from the derivative of Q(x\), producing: 

II [r(e a ) - r(e h )] K^(0,x 2 ,x 3 )h^dx 2 dx 3 . (21) 

When this is summed over a to obtain the total perturbation 
integral, however, J2 a K^HO, x 2 , x 3 ) = L(§,x 2 ,x 3 ) = 1 
by construction (since L — 1 on the interface). Thus, 
we obtain the surface integral of eq. as desired, where 
h^dx 2 dx 3 = h • dA. 

The analysis of the A/x( Q ' term proceeds identi- 
cally, although here the continuous field components are 
(Bi, H 2 H 3 ), but in this case it yields zero if fj, a = /j, b (as in 
the common case of non-magnetic materials where /1 is iden- 
tically 1). 



V. NUMERICAL VALIDATION 

To check the correctnessf of the perturbative analysis 
above, we performed the following numerical computation. 
We solve the full- vector Maxwell eigenproblem numerically, 
for inhomogeneous anisotropic dielectric structures, by iter- 
ative Rayleigh-quotient minimization in a planewave basis, 
using a freely available software package 0. Given an ar- 
bitrary structure, we can then evaluate the derivative of the 
eigenfrequency for a shifting interface, both by the perturba- 
tion eq. ([3]) and by numerical differentiation of the eigenfre- 
quencies (here, differentiating a cubic-spline interpolation). 

In particular, we considered a two-dimensional photonic 
crystal ifTBI consisting of a square lattice (lattice constant a) 
of 0.4a x 0.2a dielectric blocks of a material e a surrounded 
by e b , with Gaussian bumps on one side (inset of Fig. [3V 
Here, e a and e h are chosen to be random symmetric positive- 
definite matrices with eigenvalues ranging from 2 to 12 for 
e a and from 1 to 5 for e b . On the right side of each block 
(along one of the 0.4a edges) is a Gaussian bump of height 
h(y) — he~ v l 2w , with a width w = 0.1a and amplitude 
h (where h < denotes an indentation). We then com- 
puted the lowest eigenvalue u>(A) and eigenfields E for a set 
of h values h/a e [—0.17, +0.17] , at a Bloch wavevector 



6 




numerical derivative 

O perturbation theory 




Figure 3: Numerical validation of perturbation- theory formula, ap- 
plied to compute the derivative dw/dh for a Gaussian "bump" of 
height ft on a square lattice (period a) of anisotropic-e rectangles 
(inset) with an eigenfrequency u (corresponding to A ~ 3a). Pos- 
itive/negative h indicate bumps/indentations (see lower-right/left in- 
sets for h = ±0.15a), respectively. Solid lines are numerical dif- 
ferentiation of the eigenfrequency, and dots are from perturbation 
theory. The different lines correspond to different random dielectric 
tensors e a and e b . 

k = (0,0, 0.5) • 2n /a (leading to modes with a vacuum wave- 
length A ~ 3a). Given this data, we then compared the deriva- 
tive du/dh as computed by the perturbation equation ^ com- 
pared to the derivative of a cubic-spline fit of the frequency 
data. This was repeated for six different random e a and e b . 
The results, shown in Fig. [3] demonstrate that the perturbation 
formula indeed predicts the exact slope h as expected (with 
tiny discrepancies, due to the finite resolution, too small to 
see on this graph). 



VI. APPLICATION TO SUB-PIXEL SMOOTHING 

In any numerical method involving the solution of the full- 
vector Maxwell's equations on a discrete grid or its equiva- 
lent, such as the planewave method above O or the finite- 
difference time-domain (FDTD) method |26|, discontinuities 
in the dielectric function e (and the corresponding field dis- 
continuities) generally degrade the accuracy of the method, 
typically reducing it to only linear convergence with resolu- 
tion (6] [7). Unfortunately, piecewise-continuous e is the most 
common experimental situation, so a technique to improve the 
accuracy (without switching to an entirely different compu- 
tational method) is desirable. One simple approach that has 
been proposed by several authors is to smooth the dielectric 
function, or equivalently to set the e of each "pixel" to be some 
average of e within the pixel, rather than merely sampling e 
in a "staircase" fashion (El El |T3] [27l |28l E21 (30l [3T1 . Unfortu- 
nately, this smoothing itself changes the structure, and there- 
fore introduces errors. We analyzed this situation in a recent 



paper for the FDTD method 0, and showed that the problem 
is closely related to perturbation theory: one desires a smooth- 
ing of e that has zero first-order effect, to minimize the error 
introduced by smoothing and so that the underlying second- 
order accuracy can potentially be preserved. At an interface 
between two isotropic dielectric materials, the first-order per- 
turbation is given by eq. Q, and this leads to an anisotropic 
smoothing: one averages e^ 1 for field components perpen- 
dicular to the interface, and averages e for field components 
parallel to the interface, a result that had previously been pro- 
posed heuristically by several authors l6l [T3l |29l . 

In this section, we generalize that result to interfaces be- 
tween anisotropic materials, and illustrate numerically that it 
leads to both dramatic improvements in the absolute magni- 
tude and the convergence rate of the discretization error. In 
the anisotropic-interface case, a heuristic subpixel smoothing 
scheme was previously proposed |6|, but we now show that 
this method was suboptimal: although it is better than other 
smoothing schemes, it does not set the first-order perturba- 
tion to zero and therefore does not minimize the error or per- 
mit the possibility of second-order accuracy. Specifically, as 
discussed more explicitly below, a second-order smoothing is 
obtained by averaging r(e) and then inverting r(e) to ob- 
tain the smoothed "effective" dielectric tensor. Because this 
scheme is analytically guaranteed to eliminate the first-order 
error otherwise introduced by smoothing, we expect it to gen- 
erally lead to the smallest numerical error compared to com- 
peting smoothing schemes, and there is the hope that the over- 
all convergence rate may be quadratic with resolution. 

First, let us analyze how perturbation theory leads to a 
smoothing scheme. Suppose that we smooth the underlying 
dielectric tensor e(x) into some locally averaged tensor e(x), 
by some method to be determined below. This involves a 
change Ae = e — e, which is likely to be large near points 
where e is discontinuous (and, conversely, is zero well in- 
side regions where e is constant). In particular, suppose that 
we employ a smoothing radius (defined more precisely be- 
low) proportional to the spatial resolution Ax of our numer- 
ical method, so that Ae is zero [or at most 0(Ax 2 )] except 
within a distance ~ Ax of discontinuous interfaces. To eval- 
uate the effect of this large perturbation near an interface, we 
must employ an equivalent reformulation of eq. (|3j: 

Alo ~ J F* • Ar-Fd 3 x, 

where At — -r(e) — t(s). It is sufficient to look at the pertur- 
bation in ui, since (as we remarked in Sec.[II]i the same integral 
appears in the perturbation theory for many other quantities 
(such as scattered power, etc.). If we let xi denote the (local) 
coordinate orthogonal to the boundary, then the x\ integral is 
simply proportional to ~ J At dx\ + 0(Ax 2 ) : since F is 
continuous and At = except near the interface, we can pull 
F out of the xi integral to lowest order. That means, in or- 
der to make the first-order perturbation zero for all fields F, 
it is sufficient to have § At dx\ = 0. This is achieved by 
averaging r as follows. 

The most straightforward interpretation of "smoothing" 
would be to convolve e with some localized kernel s(x), 
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where J s(x) d 3 x = 1 and s(x) = for |x| greater than some 
smoothing radius (the support radius) proportional to the res- 
olution ~ Ax. That is, e(x) = e * s = J e(y) s(x — y) d 3 y. 
For example, the simplest subpixel smoothing, simply com- 
puting the average of e over each pixel, corresponds to s = 1 
inside a pixel at the origin and s = elsewhere. However, this 
will not lead to the desired J At = to obtain second order 
accuracy. Instead, we employ: 

e(x) - r- l [r{e) * s] = r- 1 || t [e(y)] S (x - y) d 3 yj ■ 

(22) 

where t _1 is the inverse of the r(e) mapping, given by: 

1 T12 T 13 



, Til Til TH 



Til 
.131 



<32 



Til 
T31 T 12 



-"33 



Til 
T31T13 



(23) 

The reason why eq. ( |22] i works, regardless of the smoothing 
kernel s(x), is that 



Ardh 



dh 



■ [e(y)] s(x - y) d 3 y - e(x) 



J d 3 yr[£(y)]|| s(x-y)d 3 



0. 



(24) 



This guarantees that the integral of At is zero over all space, 
but above we required what appears to be a stronger condition, 
that the local, interface-perpendicular integral J At dx\ be 
zero (at least to first order). However, in a small region where 
the interface is locally flat (to first order in the smoothing ra- 
dius), Ax must be a function of x\ only by translational sym- 
metry, and therefore ( 24 1 implies that J At dx\ = by itself. 
Although the above convolution formulas may look compli- 
cated, for the simplest smoothing kernel s(x) the procedure is 
quite simple: in each pixel, average t(s) in the pixel and then 
apply t _1 to the result. (This is not any more difficult to apply 
than the procedure implemented in Ref. [6 1, for example.) 

Strictly speaking, the use of this smoothing does not guar- 
antee second-order accuracy, even if the underlying numeri- 
cal method is nominally second-order accurate or better. For 
one thing, although we have canceled the first-order error due 
to smoothing, it may be that the next-order correction is not 
second-order. Precisely this situation occurs if one has a struc- 
ture with sharp dielectric corners, edges, or cusps, as dis- 
cussed in Ref. [1|: in this case, smoothing leads to a con- 
vergence rate between first order (what would be obtained 
with no smoothing) and second order, with the exponent de- 
termined by the nature of the field singularity that occurs at 
the corner. 



A. Numerical smoothing validation 

As a simple illustration of the efficacy of the subpixel 
smoothing we propose in eq. (22i, let us consider a two- 
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dimensional example problem: a square lattice (period a) of 



resolution (pixels/a) 

Figure 4: Relative error Alo/uj for an eigenmode calculation with 
a square lattice (period a) of 2d anisotropic ellipses (green inset), 
versus spatial resolution, for a variety of sub-pixel smoothing tech- 
niques. Straight lines for perfect linear (black dashed) and perfect 
quadratic (black solid) convergence are shown for reference. 



ellipses made of e a surrounded by e b , where we will find the 
lowest-w Bloch eigenmode. As above, we choose the dielec- 
tric tensors to be random positive-definite symmetric matrices 
with random eigenvalues in [2, 12] for e a and in [1, 5] for e b , 
and the ellipses are oriented at an arbitrary angle, at an ar- 
bitrary Bloch wavevector ka/27r = (0.1,0.2,0.3), to avoid 
fortuitous symmetry effects. (The vacuum wavelength A cor- 
responding to the eigenfrequency u> is A = 5.03a.) For each 
resolution Ax, we assign an e to each pixel by computing t _1 
of the average of r(e) within that pixel. Then, we compute 
the relative error Aw /ui (compared to a calculation at a much 
higher resolution) as a function of resolution. For compar- 
ison, we also consider four other smoothing techniques: no 
smoothing, averaging e in each pixel 1281 . averaging e _1 in 
each pixel, and a heuristic anisotropic averaging proposed by 
Ref. j6) in analogy to the scalar case. The results are shown in 
Fig.|4] based on the same planewave method as above Q, and 
show that the new smoothing technique clearly leads to the 
lowest errors Alu/lu. Also, whereas the other methods yield 
clearly first-order convergence, the new method seems to ex- 
hibit roughly second-order convergence. The no-smoothing 
case has extremely erratic errors, as is typical for stair-casing 
phenomena. 

In Fig. [5] we also show results from a similar calculation in 
three dimensions. Here, we look at the lowest eigenmode of a 
cubic lattice (period a) of 3d ellipsoids (oriented at a random 
angle) made of e a surrounded by e b , both random positive- 
definite symmetric matrices as above. The frequency u, at an 
arbitrarily chosen wavevector ka/27r = (0.4, 0.3, 0.1), corre- 
sponds to a vacuum wavelength A = 3.14a. Again, the new 
method almost always has the lowest error by a wide margin, 
especially if the unpredictable dips of the no-smoothing case 
are excluded, and is the only one to exhibit (apparently) better 
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Figure 5: Relative error Alo/uj for an eigenmode calculation with cu- 
bic lattice (period a) of 3d anisotropic ellipsoids (green inset), versus 
spatial resolution, for a variety of sub-pixel smoothing techniques. 
Straight lines for perfect linear (black dashed) and perfect quadratic 
(black solid) convergence are shown for reference. 



than linear convergence. 

Our previous heuristic proposal from Ref. [6], while bet- 
ter than the other smoothing schemes (and less erratic than no 
smoothing), is clearly inferior to the new method. Previously, 
we had observed what seemed to have been quadratic conver- 
gence from the heuristic scheme 0, but this result seems to 
have been fortuitous — as we demonstrated recently, even non- 
second-order schemes can sometimes appear to have second- 
order convergence over some range of resolutions for a partic- 
ular geometry 0. The key distinction of the new scheme, that 
lends us greater confidence in it than one or two examples can 
convey, is that it is no longer heuristic. The new smoothing 
scheme is based on a clear analytical criterion — setting the 
first-order perturbative effect of the smoothing to zero — that 
explains why it should be an accurate choice in a wide variety 
of circumstances. 



VII. CONCLUDING REMARKS 

We have shown how to correctly treat lowest-order per- 
turbations to a boundary between two anisotropic materials, 
a problem for which previous approaches had been stymied 
by the complicated discontinuous boundary conditions on the 
electric field. This result immediately led to an improved sub- 
pixel smoothing scheme for discretized numerical methods — 
we demonstrated it for a planewave method, but we expect 
that it will similarly be applicable to other methods, e.g. 
FDTD |5|. The same result can also be applied to con- 
structing an effective-medium theory for subwavelength mul- 
tilayer films of anisotropic materials. Moreover, in the pro- 
cess of deriving our perturbative result, we developed a local 
coordinate-transform approach that may be useful in treating 



many other types of interface perturbations, because it circum- 
vents the difficulty of shifting discontinuities without requir- 
ing one to construct an explicit coordinate transformation. 



Acknowledgments 

This work was supported in part by Dr. Dennis Healy of 
DARPA MTO, under award N000 14-05 -1-0700 administered 
by the Office of Naval Research. We are also grateful to 
I. Singer at MIT for helpful discussions. 



Appendix: Coordinate transformation of Maxwell's equations 



As discussed in Sec. Ill A above, any differentiable coordi- 
nate transformation of Maxwell's equations can be recast as 
merely a transformation of e and /i, with the same solutions 
E and H only multiplied by a matrix in addition to the coor- 
dinate change lfl4l . This result has been exploited by Pendry 
et al. to obtain a number of beautiful analytical results from 
cylindrical superlenses [3| to "invisibility" cloaks |4|. It (and 
related ideas) can be used to derive coupled-mode expressions 
for bending loss in optical waveguides lfT71[T9l . A similar re- 
sult has also been employed to design perfectly-matched lay- 
ers (PML), via a complex coordinate stretching, to truncate 
numerical grids [32]. It is likely that there are many other 
applications, as well as equivalent derivations, that we are not 
aware of. Here, we review the proof in a compact form, gener- 
alized to arbitrary anisotropic media. (Most previous deriva- 
tions seem to have been for isotropic media in at least one 
coordinate frame lfl4l . or for coordinate transformations with 
purely diagonal Jacobians J' where Ja depends only on x. L 
l32l . or for constant affine coordinate transforms [33].) 

We begin with the usual Maxwell's equations for Euclidean 
space (in natural units): 



V x H 

V x E 

V • (e • E) 

V • (/. • H) 



<9E 



P 
0, 



(25) 

(26) 

(27) 
(28) 



where J and p are the usual free current and charge densities, 
respectively. We will proceed in index notation, employing 
the Einstein convention whereby repeated indices are summed 



over. Ampere's Law, eq. (25 1, is now expressed: 



daHbCabc — £ 



cd ' 



dEd 
dt 



J, 



(29) 



where e a bc is the usual Levi-Civita permutation tensor and 
d a = d/dx a . Under a coordinate change x 1— ► x', if we let 

Jab — jr^r be the (non-singular) Jacobian matrix associated 

ax h 

with the coordinate transform (which may be a function of x), 
we have 



d a — Jbad'h- 



(30) 



9 



Furthermore, as in eqs. (|5]46]), let 



E a — JbaE' b , 
Ha = JbaH' b . 



Hence, eq. (|29|) becomes 



oe; 



•Jiad'iJjbEtjCabc — £cdJld—Qj- + J, 



(31) 
(32) 



(33) 



Here, the 3 la d'i — d a derivative falls on both the Jj b and 
Hj terms, but we can eliminate the former thanks to the e a b c : 
daJjb^abc = because d a J jb = d b Jj a . Then, again multi- 
plying both sides by the Jacobian J kc , we obtain 



JkcJjbJiad'iH'jeahc — Jkc^cdJld—^T + JkcJ c (34) 



d_E[ 
dt 



Noting that J ia JjbJkc^abc = e ijk det J by definition of the 
determinant, we finally have 



1 



det J 

or, back in vector notation, 



JkcZcdJld- 



dt det J 



(35) 



where J' = C ■ 3/ det J' . Thus, we see that we can interpret 
Ampere's Law in arbitrary coordinates as the usual equation 



in Euclidean coordinates, as long as we replace the materials 
etc. by eqs. |5]-[7]i. By an identical argument, we obtain 



V'xE' = - 



JiiJ T dW 



det J 



dt 



(37) 



which yields the corresponding transformation |8]l for fi. 

The transformation of the remaining divergence equations 
into equivalent forms in the new coordinates is also straight- 
forward. Gauss's Law, eq. ( |27| ), becomes 

p = d a e a bEb = Jiad , i e a bJjbE'j = Jiad[(det J~)J al ^ e'^E'^ 



= (det + (d a J~f} det J)^ 



-13 3 
(det JW^E'j, 



(38) 



which gives V' ■ (e' ■ E') = p' for p' = pj det Similarly 



for eq. (28 1. Here, we have used the fact that 



^cl^J^ det J7* — da^anmtkijiJiniJjrn/'^ 0, (39) 



from the cofactor formula for the matrix inverse, and recalling 
that daJjb^abc = from above. In particular, note that p = 
p' = and J = <^=^ J' = 0, so a non-singular 
coordinate transformation preserves the absence (or presence) 
of sources. 
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